home *** CD-ROM | disk | FTP | other *** search
/ Gekkan Dennou Club 147 / Gekkan Dennou Club - 2000.8 Vol. 147 (Japan).7z / Gekkan Dennou Club - 2000.8 Vol. 147 (Japan) (Track 1).bin / html / docs / complex / julia3.txt < prev    next >
Encoding:
Text File  |  2000-07-11  |  5.7 KB  |  254 lines
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <setjmp.h>
  4. #include <iocslib.h>
  5. #include <doslib.h>
  6. #include <math.h>
  7. #include <time.h>
  8. #include "complex.h"
  9.  
  10. #define USAGE "\
  11. 3次方程式のジュリア集合を描く\n\
  12. 使用法: julia3 [reMin imMax scLen [scSize [mxCnt [rot_f_limit [Re_A Im_A [Re_B Im_B [Re_C Im_C]]]]]]]\n\
  13.   Newton法の初期値について reMin=実部最小[-2.0], imMax=虚部最大[2.0], scLen=辺の長さ[4.0]\n\
  14.   scSize=画面サイズ(0…256×256ドット, 1…512×512ドット)[0]\n\
  15.   mxCnt=回数制限(1～65535)[256]\n\
  16.   rot_f_limit=収束半径(収束したと判断する|f(z)|の上限)[0.01]\n\
  17.   Re_A=2次の係数の実部[0.0], Im_A=2次の係数の虚部[0.0]\n\
  18.   Re_B=1次の係数の実部[0.0], Im_B=1次の係数の虚部[0.0]\n\
  19.   Re_C=定数項の実部[-1.0], Im_C=定数項の虚部[0.0]\n\
  20.     f(z)=z^3+A*z^2+B*z+C\n\
  21. "
  22.  
  23. /* 画面サイズのデフォルト */
  24. #define SC_SIZE 0
  25.  
  26. /* 計算範囲のデフォルト */
  27. #define RE_MIN (-2.0)
  28. #define IM_MAX (2.0)
  29. #define SC_LEN (4.0)
  30.  
  31. /* 最大ループ回数のデフォルト */
  32. #define MX_CNT 256
  33.  
  34. /* 収束判定の制度 */
  35. #define ROT_F_LIMIT (0.01)
  36.  
  37. /* キーチェック */
  38. #define keysns() \
  39. ({ \
  40.   int _k_ = 0; \
  41.   while (B_KEYSNS()) _k_ = B_KEYINP() & 0xff ? : _k_; \
  42.   _k_; \
  43. })
  44.  
  45. double reMin = RE_MIN, imMax= IM_MAX;
  46. double scLen = SC_LEN;
  47. int scSize = SC_SIZE;
  48. int mxCnt = MX_CNT;
  49. double rot_f_limit = ROT_F_LIMIT;
  50. double rot_f_limit2;
  51.  
  52. unsigned short *paTable;  /* パレットテーブル */
  53.  
  54. void make_palet_table(void);  /* パレットテーブルを作る */
  55. void init_screen(void);  /* 画面の初期化 */
  56. void tini_screen(void);  /* 画面の後始末 */
  57.  
  58. /* 係数と定数項のデフォルト */
  59. #define Re_A (0.0)
  60. #define Im_A (0.0)
  61. #define Re_B (0.0)
  62. #define Im_B (0.0)
  63. #define Re_C (-1.0)
  64. #define Im_C (0.0)
  65.  
  66. /* 係数と定数項 */
  67. complex A, B, C;
  68.  
  69. /* f(z),df_f(z)を変更することで別のジュリア集合を描くことができる */
  70. /* f(z)=z^3+A*z^2+B*z+C, df_f(z)=3*z^2+2*A*z+B */
  71. /* マクロ版と関数版のどちらか速い方を選択する */
  72. #if 0
  73. #define f(z) cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul((z), B)), C)
  74. #define df_f(z) cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul((z), A), 2.0)), B)
  75. #else
  76. complex f(complex z)
  77. {
  78.   return cadd(cadd(cadd(cpow3(z), cmul(cpow2(z), A)), cmul(z, B)), C);
  79. }
  80. complex df_f(complex z)
  81. {
  82.   return cadd(cadd(cremul(cpow2(z), 3.0), cremul(cmul(z, A), 2.0)), B);
  83. }
  84. #endif
  85.  
  86. /* 初期値cでf(z)に対するNewton法の漸化式の反復回数を返す */
  87. int iteration(complex c, int limit_cnt)
  88. {
  89.   complex z;
  90.   int cnt;
  91.  
  92.   if (setjmp(cjmpenv)) {
  93.     return limit_cnt;  /* 計算できなかった */
  94.   }
  95.   cnt = limit_cnt;
  96.   z = c;
  97.   while (cnt > 0 && cabs2(f(z)) > rot_f_limit2) {
  98.     /* f(z)に対するNewton法の漸化式 rec_f(z)=z-f(z)/df_f(z) */
  99.     z = csub(z, cdiv(f(z), df_f(z)));
  100.     cnt--;
  101.   }
  102.   return limit_cnt - cnt;
  103. }
  104.  
  105. /* 描画ルーチンの本体 */
  106. void draw_body()
  107. {
  108.   unsigned short *a;
  109.   complex c;
  110.   double step;
  111.   unsigned short scCnt;
  112.   unsigned short reCnt;
  113.   unsigned short imCnt;
  114.   complex c0;
  115.  
  116.   scCnt = (scSize ? 512 : 256);
  117.   step = scLen / (double)scCnt;
  118.  
  119.   /* VRAMの書き込みアドレス */
  120.   a = (unsigned short *)0xc00000;
  121.  
  122.   c0 = tocomplex(reMin, imMax);
  123.  
  124.   /* 虚軸方向のループ */
  125.   for (imCnt = 0; imCnt < scCnt; imCnt++) {
  126.     {
  127.       unsigned short *b;
  128.       b = a;
  129.       /* これから描画するラスタを示すための点線を描く */
  130.       for (reCnt = 0; reCnt < 512; reCnt += 8) {
  131.     *b = 0xffff;
  132.     b += 8;
  133.       }
  134.     }
  135.  
  136.     c = c0;
  137.  
  138.     /* 解像度に応じて実軸方向のループ全体を分ける */
  139.     if (scSize == 0) {
  140.       /* 256×256ドット */
  141.  
  142.       /* 実軸方向のループ */
  143.       for (reCnt = 0; reCnt < scCnt; reCnt++) {
  144.     unsigned short col;
  145.  
  146.     /* 漸化式の反復回数を数えてパレットテーブルを参照する */
  147.     col = paTable[iteration(c, mxCnt)];
  148.  
  149.     /* 解像度に応じた大きさの点をVRAMに書き込む */
  150.     *a++ = col;
  151.     *a++ = col;
  152.     a[512-2] = col;
  153.     a[512-1] = col;
  154.  
  155.     /* 実軸方向に進む */
  156.     Re(c) += step;
  157.       }
  158.     } else {
  159.       /* 512×512ドット */
  160.  
  161.       /* 実軸方向のループ */
  162.       for (reCnt = 0; reCnt < scCnt; reCnt++) {
  163.     /* 漸化式の反復回数を数えてパレットテーブルを参照する */
  164.     /* 解像度に応じた大きさの点をVRAMに書き込む */
  165.     *a++ = paTable[iteration(c, mxCnt)];
  166.  
  167.     /* 実軸方向に進む */
  168.     Re(c) += step;
  169.       }
  170.     }
  171.  
  172.     /* 解像度に応じてVRAMの書き込みアドレスを補正する */
  173.     if (scSize == 0) {
  174.       a += 512;
  175.     }
  176.  
  177.     /* 虚軸方向に進む */
  178.     Im(c0) -= step;
  179.  
  180.     /* キーが押されたら中止 */
  181.     if (keysns()) {
  182.       break;
  183.     }
  184.   }
  185. }
  186.  
  187. void main(int argc, char *argv[])
  188. {
  189.   clock_t tm0, tm1;
  190.  
  191.   /* 画面の初期化 */
  192.   init_screen();
  193.  
  194.   /* 係数と定数項のデフォルトを設定 */
  195.   A = tocomplex(Re_A, Im_A);
  196.   B = tocomplex(Re_B, Im_B);
  197.   C = tocomplex(Re_C, Im_C);
  198.  
  199.   /* コマンドラインのパラメータを取得 */
  200.   switch (argc) {
  201.   case 13:
  202.     sscanf(argv[12], "%lf", &Im(C));
  203.     sscanf(argv[11], "%lf", &Re(C));
  204.   case 11:
  205.     sscanf(argv[10], "%lf", &Im(B));
  206.     sscanf(argv[9], "%lf", &Re(B));
  207.   case 9:
  208.     sscanf(argv[8], "%lf", &Im(A));
  209.     sscanf(argv[7], "%lf", &Re(A));
  210.   case 7:
  211.     sscanf(argv[6], "%lf", &rot_f_limit);
  212.   case 6:
  213.     sscanf(argv[5], "%d", &mxCnt);
  214.   case 5:
  215.     sscanf(argv[4], "%d", &scSize);
  216.   case 4:
  217.     sscanf(argv[3], "%lf", &scLen);
  218.     sscanf(argv[2], "%lf", &imMax);
  219.     sscanf(argv[1], "%lf", &reMin);
  220.   case 1:
  221.     break;
  222.   default:
  223.     printf(USAGE);
  224.     return;
  225.   }
  226.  
  227.   /* 収束半径を2乗しておく */
  228.   rot_f_limit2 = rot_f_limit * rot_f_limit;
  229.  
  230.   /* パレットテーブルを作る */
  231.   make_palet_table();
  232.  
  233.   /* VRAMに直接書き込むためにスーパーバイザモードにする */
  234.   B_SUPER(0);
  235.  
  236.   /* 計測開始 */
  237.   {
  238.     clock_t t = clock();
  239.     while ((tm0 = clock()) == t);
  240.   }
  241.  
  242.   /* 描画ルーチンの本体 */
  243.   draw_body();
  244.  
  245.   /* 計測終了 */
  246.   tm1 = clock();
  247.  
  248.   /* 画面の後始末 */
  249.   tini_screen();
  250.  
  251.   /* 所要時間を表示 */
  252.   printf("所要時間は %.8g 秒です．\n", ((double)(tm1 - tm0)) / CLK_TCK);
  253. }
  254.